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Introduction 


A  method  for  generating  a  normal  random  variable  in  terms  of  uniform 
random  variables  will  be  described.  The  method  is  based  on  representing 
a  density  function  as  a  mixture  of  simpler  densities,  as  outlined  in  [ 1 ] . 
It  is  quite  fast  and  requires  little  storage  (60  constants).  It  is  not 
quite  as  fast  as  the  method  of  [2],  but  it  is  simpler,  with  less  chance 
for  prospective  users  being  set  adrift  in  a  sea  of  details. 

We  will  fashion  the  procedure  with  use  on  a  computer  with  numbers 
having  35  binary  digits  (bits)  in  mind.  The  details  may  be  modified  for 
numbers  of  different  lengths  or  for  decimal  computers. 

The  normal  random  variable  x  is  generated  in  terms  of  independent 
uniform  random  variables  u^jU^,...  .  Roughly,  the  procedure  is  this: 
we  put  x  =  ^u  with  frequency  35%;  x  =  £  +  ^u,  24%;  x  =  1  +  -^u,  13%; 
x  =  3/?  +  ^u,  5%;  and  x  =  2  +  2U,  with  frequency  about  1%.  Thus 
with  probability  about  .8  we  put  x  =  a  +  ^u,  the  first  9  bits  of  u^ 
being  used  to  choose  a  from  among  2,  l/2,  1,  3/2,  2  and  the 
remaining  26  bits  of  u^  used  to  form  ^u. 

Of  the  remaining  20%  we  use  the  modified  rejection  technique  given 
in  [3]  with  frequency  19%  and  the  other  1%  is  devoted  to  the  tail.  The 
modified  rejection  technique  is  applied  to  one  of  regions  6-10  above 
the  rectangles  in  Fig.  1,  while  rectangles  1-5  provide  the  representations 
x  =  0  +  ^u,  £  +  ^u,  etc.  We  use  the  polar  representation  of  a  random 
normal  point  to  dispose  of  the  tail. 
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The  procedure  has,  in  theory,  unlimited  accuracy;  in  practice,  its 
accuracy  is  limited  by  the  capabilities  of  the  machine  in  single 
precision  floating  point. 

2.  The  Rectangles 

Let  f  be  the  absolute  normal  density, 

f(x)  =  -  e  2  ,  x  >  0. 

\f  2tt 

We  generate  x  with  density  f,  attaching  a  random  +  at  some  convenient 
point  in  the  program.  We  write  f  as  a  mixture  of  11  densities : 

11 

f(x)  =  Z  p.g.  (x). 
i=l 

Fig.  1  shows  how  the  terms  in  the  mixture  are  stacked  to  form  f. 

The  5  rectangles  provide  the  fast  part  of  the  program;  g^,...,g^ 
being  rectangles  with  base  -g-,  an  x  with  such  a  density  may  be  found 
by  putting  x  =  a  +  -£u.  For  example,  region  1  of  Fig.  1  provides  the 

representation  x  =  0  +  -gu,  and  this  should  be  used  with  frequency  equal 

to  the  area  of  region  1,  p^  =  .352065326764299.  But  it  is  wasteful  to 
test  u^  <  p^  in  its  full  form,  as  this  will  require  all  35  bits  of  u^; 

much  better  is  to  use  only,  say,  9  bits  of  u^  to  test  u^  <  p^  where 

is  the  nearest  we  can  get  to  p^  using  9  bits,  p^  =  180/512  =  .3515625, 
then  we  may  use  the  last  26  bits  of  u^  to  form  -^u  and  put  x  =  Ju. 

Later,  with  probability  p^  -  p^,  we  will  put  x  =  ^u^,  using  a 
completely  new  uniform  number  u^.  We  then  will  get  by  most  of  the  time 
with  a  single  uniform  number  u^.  Each  of  the  probabilities 
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is  broken  up  in  this  way;  the  numerical  values  appear  in  the  steps  of 
the  program  and  in  the  summary  of  constants. 


3*  The  "Nearly  Linear11  Method 


We  illustrate  the  nearly  linear  method  with  region  9  of  Fig.  1. 

It  represents  the  term  Pggg  in  the  mixture,  so  we  must,  with 
probability  pg  =  .034123172128170,  generate  a  random  variable  with 
density  gg.  This  density  is  drawn  in  Fig.  2.  Our  method  is  to  enclose 
the  area  under  g  in  the  trapezoid  ADEF.  The  slope  of  CG  and  DE 
is  -8.  We  choose  a  point  (z,y)  uniformly  from  the  trapezoid  and  set 
x  =  z  if  y  <  g(z) ,  but  choose  a  new  point  (z,y)  if  not.  We  test  to 
see  if  the  point  (z,y)  is  in  triangle  ACG  (y  <  const.  -8z),  hoping  to 
avoid  testing  y  <  g(z).  If  that  test  is  necessary,  we  take  advantage 
of  the  exponential  content  of  g,  writing 


g(z)  =  [e-U^2-9/i)  _  e-7/Sh 

Pg  P9  V  2tt 

Then  the  test  y  <  g(z)  hinges  on  whether  this  infinite  sum  is  positive 
or  negative: 


ry 


s  +  t 


with  r  =  pg  sIV/2  e^^,  s  =  1  -  ,  t  =  -g{z  +  3/2)  (z  -  3/2), 

0  <  t  <  7/8.  This  is  an  alternating  series  that  converges  rather 
quickly,  and  its  sum  is  negative  if  and  only  if  two  successive  partial 
sums  are  negative;  positive  if  and  only  if  two  successive  partial  sums 
are  positive.  We  should  average  only  a  few  multiplications  in  order  to 


decide  y  <  g(z)  with  "unlimited"  accuracy.  Another  advantage  is  that 
the  structure  of  the  test  for  gg  is  the  same  as  that  for  the  other  g's, 
only  r,  s,  and  t  changing. 

The  point  (z,y)  in  the  trapezoid  of  Fig.  2  may  be  produced  quite 
fast,  putting 

z  =  minCu^jU^)  ■+*  *-5» 

y  =  d  +  *fmax  (up,^) 

(choosing  from  triangle  BDE)  with  probability  and  putting 

2  =  £u2, 

y  =  du^ 

(choosing  from  ABEF)  with  probability  ;r—  •  If  d  is  small  we  will 
avoid  conventional  multiplication  most  of  the  time,  and  by  properly 
choosing  d  we  will  need  only  a  few  bits  of  u^  to  choose  between  the 
two  methods  for  generating  (z,y),  so  that  the  remaining  bits  of  u^  may 
be  used  to  form  Up  and  u^« 

For  example,  pg  =  .034123172128170.  We  let  pg  =  17/512  =  .033203125, 
then  generate  x  with  density  gg  with  probability  pg.  In  the  course 
of  the  program,  the  test  for  this  contingency  is  <  u^  <  • 

If  u^  is  in  this  range,  then  by  making  ,  we  may  test 

U1  ^  ^^512^  ’  keePinS  the  26  bits  of  u^  available  for  forming 

u?  and  ug.  This  makes  d  =  3/7  =  .42857...  whereas  the  smallest 
possible  value  for  d  is  d  =  .427...  .  The  loss  is  small,  and  the  gain 
from  not  having  to  generate  new  Up  and  ug  is  great.  We  must  generate  x 
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from  density  gQ  with  probability  pQ  -  pn  later  on,  using  new  u„ 
y  y  y  <Z 

and  u^>  but  the  average  time  loss  for  this  compensating  step  is 
negligible. 

4.  The  Tall 

We  need  x  conditioned  by  jx|  >  2.5.  Consider  the  plane  of  points 
(x-^x^)  with  x^  and  x^  independent  standard  normal. 


We  want  a  point  (x^jX^)  outside  the  square,  then  one  of  the  coordinates 
will  serve  as  x.  We  use  the  polar  representation  (p,0)  to  choose  points 
outside  the  circle  until  we  get  one  outside  the  square,  then  change  to 
rectangular  coordinates  (x^,x^). 

The  exterior  of  the  circle  has  measure  .044,  the  exterior  of  the 
square,  .025.  The  ratio  of  the  two  is  about  .57  so  that  57%  of  the  time 
a  point  lying  outside  the  circle  will  lie  outside  the  square.  If  (x^,x?) 
lies  outside  the  square,  only  once  in  150  times  will  both  |x^|  and  | x^ I 
be  >  2.5,  so  it  probably  isn't  worth  the  trouble  to  store  the  extra  one 
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for  future  use,  as  this  would  mean  that  each  time  the  tail  was  encountered 
in  the  program,  the  request  for  a  possible  stored  value  of  x  would  almost 
always  be  denied . 

It  is  elementary  to  verify  that  putting 


with  w  exponential  and  v-^v^  uniform  on  (-1,1)  and  conditioned  by 
2  2 

v^  +  v^  <  1,  will  provide  a  point  (x^,x^)  outside  the  circle  with  the 
appropriate  distribution. 

5 .  Outline  of  a  Program 

We  give  a  series  of  steps  outlining  a  program  based  on  the  above 

procedure,  assuming  we  have  a  subroutine  which,  when  called  for,  will 

generate  a  new  uniform  number  u  in  terms  of  the  current  u;  for 

k  35 

example,  by  putting:  new  u  =  2  u  +u  +  1  (mod  2JJ) ,  as  suggested  by 
Rotenberg  [5 ] . 
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i  an 

1.  If  0  <  u^  <  put  x  =  gu,  with  -gu  formed  from  bits  10-35  of  u^. 

2.  If  P11^  x  =  -g-  +  -£u,  with  formed  from  bits  10-35  of  u^. 

3.  If  put  x  =  1  +  gu,  with  formed  from  bits  10-35  of  u^. 

4.  If  ^2"  ^  ^  $12  put  x  =  1.5  +  gu*,  with  -gu  formed  from  bits  10-35  of  u^. 

5.  If  |2|  <  u^  <  put  x  =  2.0  +  -gu,  with  -gu  formed  from  bits  10-35  of  u-^. 

6*  If  512-ul<5lf  ^  j  =  7»  if  512^ul<5lf  put  j  =  8> 

if  fit  ^  U1  <  512  17111  -j  =  9»  if  512  ^  ui  <  5?f  ^  i  =  6>  and 

if  <  u^  <  pit  j  =  10.  Test  u^  <  q^,  if  yes,  form 

z  =  a.  +  jr  min(u0,u_)  and  y  =  d.  +  4max(u.,u„);  if  no,  form 
•  J  <■  3  J  2  3 


z  =  +  ^u2,  y  =  d^u^.  In  either  case,  u^  is  formed  from 

bits  10-22  and  u^  from  bits  23-35  of  u^.  [K]  If  y  <  e.  -  8z, 

put  x  =  zj  if  not,  let  t  =  j?{z  +  a.)(z  -  a.)  and  form  terms  of 

J  <J 


the  sequence 


W  r/-sj+t’  r/-»j+t-2T> 


r  .y-s  ,+t 
r  J 


until  two  successive  terms  have  the  same  sign.  If  two  successive 
terms  are  negative,  put  x  =  z,  if  two  successive  terms  are  posi¬ 
tive,  go  to  H. 

7.  If  <  u^  <  .980619788956091  choose  i  from  1,2,..., 5  with 
probability  -  p^ ,  generate  a  uniform  number  u^  and  put 
x  -  -g{i-l)  +  ^u2> 


\ 
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8.  If  .980619788956091  <  u^  <  .987580669348448  choose  j  from  among 

6,7,8,9,10  with  probability  p.  -  p.  and  go  to  H. 

J  J 

9.  If  .987580669348448  <  u,  <1,  form 


where  w  has  the  exponential  distribution  and  v-^v^  are  unif'orm 

o  o 

on  (-1,1),  conditioned  by  v^  +  v^  <  1.  Test  |x^|  >  2.5;  if 
yes,  put  x  =  x^.  If  no,  test  Ix^J  >  2.5;  if  yes,  put  x  =  x^; 
if  no,  generate  a  new  pair  x^,x,,  and  try  again. 


H. 


Generate  a  new  uniform  number  u.  Test  u  <  b.,  if  yes,  let 

J 

z  =  a^  +  g-  min(u2,u^) ,  y  =  d^  +  4max(u?,u^);  if  no,  let 
z  =  a^  +  gu,,,  y  =  dju3  where  u?  and  u^  are  formed  from  the 
first  and  last  parts  of  a  new  uniform  number  generated  from  u. 
Now  go  to  [K]  of  step  6. 


6.  Remarks 

With  probability*  =  .79  we  will  use  the  representation 

x  =  a  +  -gu  of  steps  1-5,  and  this  is  the  principal  reason  for  a  fast 

95  r* 

program.  However,  even  when,  with  frequency  or  about  19%,  we 

require  step  6,  the  time  loss  is  not  crushing,  for  most  of  the  time  we 
will  get  through  step  6  without  using  a  single  conventional  multipli¬ 
cation  and  will  require  only  one  uniform  number  u^. 
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It  is  in  trying  to  avoid  the  use  of  more  than  one  uniform  number 

in  step  6  that  we  make  our  only  compromise  with  accuracy,  for  we  use 

only  13  bits  to  represent  numbers  on  an  interval  of  length  jr.  A 

program  based  on  the  outline  above  can  be  described  as  being  as  accurate 

as  the  single  precision  capacity  of  the  machine,  except  that  with 

probability  j  the  value  of  x  may  differ  by  about  2~' ^  from  the 

value  that  should  be  delivered.  Any  nuts  wanting  accuracy  greater  than 

this  can  call  up  a  new  uniform  number  at  that  point,  with  a  loss  of 

15 

perhaps  15  cycles,  or  a  net  loss  of  about  —  =  3  cycles  in  the  average 
running  time  of  the  program. 

Step  7  corrects  for  our  inability  to  handle  regions  1  to  5  using 
only  9  bits,  and  step  8  does  the  same  for  regions  6  to  10.  Step  9 
provides  an  x  from  the  tail.  Note  that  in  step  9  a  signed  value  of 
x  is  returned,  whereas  in  the  other  steps  a  positive  x  is  returned 
and  a  random  +  must  be  assigned  at  some  convenient  point.  The  method 
for  doing  this  is  not  critical  in  steps  7  or  8,  but  in  the  first  6  steps, 
particularly  steps  1-5,  the  method  should  be  chosen  so  that  time  is  not 
wasted  in  this  assignment. 

Step  9  requires  an  exponential  variable  w.  The  ease  with  which 
w  can  be  produced  is  probably  more  important  than  saving  time  here, 
since  step  9  provides  only  one  x  in  a  hundred.  If  a  logarithm  sub¬ 
routine  is  available,  w  =  -  In  u  will  serve,  or  the  method  in  [4] 


may  be  used 
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SUMMARY  OF  CONSTANTS 


j 

b. 

J 

d  . 

J 

A 

PJ 

a . 
j 

6 

12/15 

1/2 

3.0380 

489/512 

15/512=. 029296875 

0 

7 

28/29 

1/14 

7.8099 

432/512 

29/512=. 056640625 

1/2 

8 

25/27 

4/25 

11.930 

458/512 

27/512=. 052734375 

1 

9 

14/17 

3/7 

15.850 

474/512 

17/512=. 033203125 

3/2 

10 

5/7 

4/5 

19.780 

497/512 

7/512=. 013671875 

2 

j 

p; 

A 

s  . 

J 

6  .030859595783727 

7  .057793845069917 

8  .054178509659306 

9  .034123172128170 

10  .015552632751237 


.001562720783727 

.001153220069917 

.001444134659306 

.000920047128170 

.001880757751237 


.038676767667587 

.082078297231197 

.111952612794322 

.131731800427316 

.144029953116673 


.117503097415405 

.312710721209028 

.464738571481000 

.583137980321492 

.675347532641650 


1 

1 

2 

3 

4 

5 


pi 


pi 


.352065326764299 

.241970724519143 

.129517595665892 

.053990966513188 

.017528300493569 


180/512=. 3515625 
123/512=. 240234375 
66/512=. 12890625 
27/512=. 052734375 
8/512=. 015625 


.000502826764299 

.001736349519143 

.000611345665892 

.001256591513188 

.001903300493569 
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